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Abstract. We describe a method for discrctizing planar C* 2 -regular domains immersed in non- 
conforming triangulations. The method consists in constructing mappings from triangles in a back- 
ground mesh to curvilinear ones that conform exactly to the immersed domain. Constructing such 
a map relies on a novel way of parameterizing the immersed boundary over a collection of nearby 
edges with its closest point projection. By interpolating the mappings to curvilinear triangles at 
select points, we recover isoparametric mappings for the immersed domain defined over the back- 
ground mesh. Indeed, interpolating the constructed mappings just at the vertices of the background 
mesh yields a fast meshing algorithm that involves only perturbing a few vertices near the boundary. 

For the discretization of a curved domain to be robust, we have to impose restrictions on the 
background mesh. Conversely, these restrictions define a family of domains that can be discretized 
with a given background mesh. We then say that the background mesh is a universal mesh for such 
a family of domains. The notion of universal meshes is particularly useful in free/moving boundary 
problems because the same background mesh can serve as the universal mesh for the evolving domain 
for time intervals that are independent of the time step. Hence it facilitates a framework for finite 
element calculations over evolving domains while using a fixed background mesh. Furthermore, 
since the evolving geometry can be approximated with any desired order, numerical solutions can be 
computed with high-order accuracy. We demonstrate these ideas with various numerical examples. 
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1. Introduction. Finite element methods commonly handle evolving domains 
in one of two ways — either the changing domain is rcmeshed at each instant /update, 
or it is immersed in a background mesh and approximated within it. We introduce a 
novel approach here that inherits the conceptual simplicity of the former and the com- 
putational efficiency of the latter. We describe a method for discretizing sufficiently 
smooth planar domains using a given background triangulation, provided some con- 
ditions are met. The method consists in constructing mappings from triangles in a 
background mesh to curvilinear ones that conform exactly to the immersed domain. 
As an example, consider simulating a problem in which rigid blades physically mix 
fluid in a closed container. As the blades rotate, the region of the container occupied 
by the fluid changes. We can now discretize the evolving fluid domain by merely per- 
turbing vertices and edges of the same background mesh. Precisely because the same 
background mesh is utilized for all positions of the propeller, we term it a universal 
mesh for the fluid. Since connectivities of triangles remain unaltered and no new 
vertices are introduced, sparse patterns of data structures involved in the problem 
can also be retained. 

The basic idea: The key idea in constructing conforming discretizations for an 
immersed domain consists of a two-step procedure to perturb edges and vertices in 
its background mesh. We identify triangles in the background mesh with at least 
one vertex inside the curved domain. Edges belonging to this collection having both 
vertices outside the domain are mapped onto the boundary with its closest point 
projection. Then we relax away from the boundary a few vertices that lie inside the 
domain and close to the boundary. For a certain class of background meshes, these 
steps enable us to construct a homeomorphism between the union of selected triangles 
in the background mesh and the the domain. In other words, this construction yields 



a conforming curvilinear discretization for the immersed domain. 

There are no conformity requirements on the background mesh; for instance, 
none of its vertices need to lie on the boundary. The resulting algorithmic advantages 
are significant, especially for problems with evolving domains. Such problems are 
ubiquitous, including ones with interaction between fluids and solids, problems with 
free boundaries and moving interfaces, domains with propagating cracks and problems 
with phase transformations. Various numerical schemes have been proposed for such 
applications, see references [3, 10, 19, 20, 27] for a representative few. The spirit 
of this article, as evidenced by the examples presented, is to immerse such evolving 
geometries in a universal mesh and update its spatial discretization as necessary. 

Numerical methods that adopt nonconforming meshes have been formulated in 
various ways. For instance, by re-triangulating elements near the boundary and with 
cut/trimmed cells [18, 24]; by treating immersed boundaries and interfaces via con- 
straints using penalty [1], Lagrange multipliers [6] and Nitsche's method [14]; or by 
enriching the space of solutions near the boundary as done in extended finite clement 
methods and discontinuous Galerkin-based methods [27, 16]. One of the challenges 
in these methods is achieving optimal accuracy with high-order interpolations. Al- 
most without exception, these methods resort to a polygonal approximation for the 
immersed domain. Such approximations suffice in low order methods, in which the 
solution is approximated by piecewise constant or affine functions [16, 22]. To con- 
struct high-order methods, it is imperative to approximate the immersed geometry 
sufficiently well over the background mesh. 

With conforming curvilinear discretizations for immersed domains, we can con- 
struct curved finite elements with optimal convergence rates given only nonconform- 
ing background meshes (§3). Alternately, by interpolating the mappings to curvi- 
linear triangles at select points, we recover isoparametric mappings defined over the 
background mesh, as discussed in §4. Rather than discretizing a curved domain ex- 
actly, isoparametric mappings provide a systematic way of constructing sufficiently 
accurate approximations in which curved edges interpolate the boundary at select 
points/nodes. Both exactly conforming and isoparametric elements enable high-order 
convergence rates (see §3.2). The former type of elements can be particularly advan- 
tageous in problems sensitive to boundary conditions. One such problem is that of a 
simply supported circular plate in bending [2]. In §4.2, we show that the accuracy of 
its numerical solution can be sensitive to how well curved boundaries are represented. 

The construction for curved domains given here are inspired by the well known 
mappings proposed in references [13, 15, 17, 26, 28]. The constructions in these articles 
assume as a point of departure, (i) a conforming mesh for the curved domain, and (ii) a 
parametrization for the curved boundary over edges that interpolate it. By admitting 
nonconforming meshes and not relying on a specific representation for the boundary, 
we relax both assumptions. Hence we generalize these known constructions to a larger 
class of meshes and amenable to general boundary representations. For example, the 
boundary can be given parametrically as a collection of splines, or implicitly as the 
zero level set of a function. We only require a way to compute the closest point 
projection sufficiently close to the boundary. 

In the special case when the mappings to curvilinear triangles are interpolated 
just at vertices in the background mesh, we get a conforming mesh for the immersed 
domain. The resulting "meshing algorithm" is very fast, since it only involves perturb- 
ing vertices. For simplicity, we present the meshing algorithm first, in §2, by defining 
affine mappings that perturb vertices of the background mesh near the immersed 
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boundary. In contrast, conforming curvilinear discretizations, discussed in §3, are 
constructed by mapping certain edges, rather than just vertices, onto the boundary. 

It is common knowledge that perturbing vertices in a background mesh can yield 
a conforming mesh for an immersed domain. The challenge, though, is making such 
algorithms robust. An algorithm specific to the case of domains immersed in rect- 
angular grids is described in [5]. The closest point projection has also been used to 
locally modify Cartesian grids near the boundary, as done in [12, 25] although the 
authors do not consider when their approach could fail. 

In §2.1, we list sufficient conditions for our meshing algorithm to be robust. We 
assume that the curved domain is C 2 -regular, and we have to restrict the class of 
background meshes for a given domain. In particular, we require that the background 
mesh be sufficiently refined and that certain angles in triangles near the boundary 
be acute. With these assumptions, we analyzed the restriction of the closest point 
projection to the edges whose vertices are projected onto the immersed boundary in 
[21]. We proved that this mapping is in fact a homeomorphism onto the boundary. 
As we discuss in §5, such a result and a possibly smaller mesh size ensures that 
moving vertices of the background mesh in the way we do in the meshing algorithm 
will not result in degenerate, inverted, or overlapping triangles, and in general, avoids 
tangled meshes. For instance, a refined background mesh of equilateral triangles is 
guaranteed to mesh a smooth domain immersed in it. We cannot however make the 
same claim with a background mesh of right angled triangles, because such a mesh 
may not satisfy the condition on angles. The conditions required for the success of 
the meshing algorithm along with possibly smaller mesh size near the boundary also 
ensure that the mappings to conforming and isoparametric curved elements are well 
defined as well. 

2. From background meshes to conforming meshes. We begin by illustrat- 
ing the steps to determine a conforming mesh for a planar curved domain f2 immersed 
in a background triangulation Th, where h denotes the mesh size. We assume that fl 
is an open set and denote f2 c = R 2 \ fi. By f2 being immersed in Th, we mean that the 
set triangulated by Th contains f2. To ensure that the resulting mesh for is valid, 
we require certain assumptions on fl and Th- They are stated in §2.1 and discussed 
in §5. 

Table 2.1: Steps in the meshing algorithm 



Step 1: Identify vertices 
in the background mesh Th 
that lie in £1 (square mark- 
ers) and in fl c (circular 
markers) respectively. Omit 
triangles with no vertices in 
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Step 2: Identify positively 
cut triangles in Tk — trian- 
gles with precisely one ver- 
tex in f2. These are shaded 
in gray. In each such trian- 
gle, check that the angle at 
its vertex in Q c closest to 9f2 
is smaller than 90°. These 
angles are labeled ■& in the 
adjacent figure. 




Step 3: A positive 
is the edge of a positively 
cut triangle joining its ver- 
tices in Ct c . These edges 
are shown in black. Map 
vertices of positive edges to 
their closest point on the 
boundary 9f2. 

Step 4: Identify vertices in 
f2 that lie close to dft, such 
as the ones shown by tri- 
angular markers. Perturb 
these vertices away from <9f2 
using the mapping p h in 
(2.1). 



Final Mesh 




The mapping for relaxing vertices away from the boundary is given by 

Ph{x)= \ X -Vh(l + ^)^) * -r <*(*)<(), (21) 
[ x otherwise, 

where <fi is the signed distance function to dtt defined below in §2.4, r/ e (0, 1) and r 
equals a few multiples of the mesh size h. We discuss the choice of v and r in §2.4 as 
well. 

The meshing algorithm in Table 2.1 is succinctly summarized as a piecewise affinc 
mapping over triangles in Th- To this end, let 7r : R 2 — > dil denote the closest point 
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projection onto d£l, and identify the collections of triangles 

7^ = {K e % : 4> > at precisely i vertices of K} (2.2) 

for i = 0, 1, 2, 3. Consider a triangle K <E T^' 1 ' 2 with vertices {u, u, w} ordered such 
that 4>{u) > (f)(v) > 4>(w). Denote the barycentric coordinates of x £ K by (A u , A„, A^) 
so that x = X u u + X v v + X w w and X u + X v + X w = 1. The algorithm in Table 2.1 
maps x M> Mk(x) defined as 

X u p h (u) + X v p h (v) + X w p h (w) itKe TH , 
X u tt(u) + X v p h (v) + X w p h (w) if K e 7^, (2.3) 
X u ir(u) + X v tt(v) + X w p h (w) if K e 7^ 2 . 

To refer to the angles checked in step 2 in Table 2.1, we introduce the terms 
■proximal vertices and conditioning angles. The proximal vertex of a positively cut 
triangle is the vertex of its positive edge closer to dft. If both vertices of the positive 
edge are equidistant from dtt, the one containing the smaller interior angle is chosen 
as the proximal vertex. If the angles are equal as well, the choice is arbitrary. The 
conditioning angle of a positively cut triangle is the interior angle at its proximal 
vertex. Hence in step 2 of Table 2.1, we require that each positively cut triangle have 
an acute conditioning angle. 

2.1. Sufficient conditions for a valid mesh. We require a few assumptions 
to guarantee that the meshing algorithm in Table 2.1 yields a valid mesh for the 
immersed domain Q. By a valid mesh we mean a triangulation of a polygon Qh made 
of triangles with diameters smaller or equal to h, such that the vertices of fi^ lie on 
dQ, and f2/j approximates fi as h \ 1 . We require that 

(a) the domain il be C 2 -regular, 

(b) that O be immersed in Th, i.e., £1 C LiKeT h ^' 

(c) the conditioning angle in each positively cut triangle in Th be strictly smaller than 
90°, and 

(d) the triangulation Th be sufficiently refined in the vicinity of dCt. 

A precise definition of C fc -regular domains is given in [21]. For our purposes, it 
suffices to note that il is C 2 -regular if the signed distance function to dil is C 2 in a 
neighborhood of dQ. Assumption (d) requires that the mesh size be small near dQ. 
By this we mean that if triangle K £ Th lies near dfl, then its diameter Kk should 
be smaller than a value that depends on the local curvature and feature size of <9f2, 
among others. Explicit estimates for the required mesh size near dfl are essential to 
automate the algorithm. However, we do not provide such estimates here. We do 
briefly mention how to construct adaptively refined background meshes satisfying the 
acute conditioning angle requirement (c) in §2.4. We discuss the rationale behind 
these assumptions later in §5. 

2.2. An illustrative example. An example that uses the meshing algorithm 
is shown in Fig. 2.1. The curved domain to be meshed is the one shown in Fig. 2.1a. 
It is C 2 -regular because its boundary is a collection of cubic splines. It is immersed 
in a background mesh of equilateral triangles. Hence the conditioning angle equals 
60° and the check in step 2 in Table 2.1 (also assumption (c) in §2.1) is trivially 
satisfied. In Fig. 2.1a, the 2382 triangles in 7^°' 1,2 are shaded in gray. Triangles 



M h K {x) = { 



x At least we need to have |f2h \ Q U f2 \ — ¥ and distancc(i9f!^, dQ) ~ > as h \ 0. 
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(a) A curved domain is immersed in a background mesh of equilateral triangles. 
The boundary of the domain shown is composed of 41 cubic splines. Its is 
therefore C 2 -regular. Since all interior angles of triangles in the background 
mesh equal 60°, the conditioning angle also equals 60°. Triangles with at 
least one vertex inside the domain (i.e., triangles in T^' 1 ' 2 ) are shaded in 
gray. These triangles are mapped to a conforming mesh for the immersed 
domain. 




(b) There are 2382 triangles in the collection T^' 1 ' 2 . Of these, 234 tri- 
angles belong to (shaded in black), 228 triangles to T? (shaded in 
dark gray) and the remaining 1920 triangles to The 756 triangles 
left unshaded in the figure are the ones that are retained unaltered 
from the background mesh. The remaining ones have their vertices 
either snapped onto the boundary or relaxed away from the boundary. 

Fig. 2.1: Example illustrating the algorithm in Table 2.1 to mesh an immersed domain 
by perturbing vertices in a background mesh. 
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(a) (b) 

Fig. 2.2: The conforming mesh determined by the meshing algorithm for the domain 
and background mesh shown in Fig. 2.1. A closer view is shown in (b). 



in this collection are mapped to a conforming mesh for the immersed domain as 
shown in Fig. 2.2. Of these triangles, 756 remain unaltered from the background 
mesh. At least one vertex in the remaining triangles is perturbed. Note that with a 
more refined background mesh, a larger fraction of the triangles in 7^°' 1,2 will remain 
unaltered from the background mesh. For instance, when the background mesh in 
Fig. 2.1 is refined once by subdivision, 5552 of the 9178 triangles in the resulting 
conforming mesh remain unaltered (i.e., remain equilateral triangles). 

In Table 2.2, we inspect the quality of triangles in the mesh in Fig. 2.2. We use 
the ratio of the circumradius to the inradius as a metric for the quality of triangles. 
The best possible value of this ratio is 2, which is attained in equilateral triangles. The 
table lists the number of triangles in the final mesh with quality in a given range of the 
metric. The minimum and maximum angles in the mesh were 20.6° and 129.6°. Table 
2.2 also reports the quality of the mesh determined by the algorithm upon refining 
the background mesh in Fig. 2.1a by subdividing each triangle into four self-similar 
ones. Extreme angles in the resulting mesh for the curved domain were 18.4° and 
139.7° respectively. 

2.3. Background meshes as Universal meshes. An important advantage of 
admitting nonconforming background meshes is in problems with evolving domains. 
For then it is possible, at least in principle, to use the same background mesh to 
triangulate a changing domain. If not for the entire duration of interest, at least for 
reasonably large changes in the immersed geometry. This motivates the notion of 
universal meshes. 

Given a triangulation Th, let T>(Th) denote the class of all domains that can be 
meshed with the algorithm in Table 2.1 using Th as a background mesh. We say 
that Th is a universal mesh for domains in T>(Th)- The utility of this concept lies 
in the fact that if {fit}* is the time evolution of a domain flo € T>(Th), then often 
{fi t : < t < T} C T>(Th) for a reasonably large time T > 0. As the domain develops 
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Fig. 2.3: Example to illustrate the notion of a universal mesh. A three blade propeller 
rotates about an axis perpendicular to its plane and passing through its center. It is 
immersed in a rehned background mesh of acute angled triangles shown in (a). Using 
this background mesh in the meshing algorithm in Table 2.1 yields a conforming mesh 
for each orientation of the propeller; a few are shown in (b). Such a background mesh 
is hence termed a universal mesh for the domain of the propeller. 
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Table 2.2: Quality of the mesh determined by the meshing algorithm for the domain in 
Fig. 2.1a. The metric used for the quality of a triangle is the ratio of its circumradius 
to the inradius. We only inspect triangles in the final mesh that have been perturbed 
with respect to the background mesh by the meshing algorithm. The remaining 
triangles remain equilateral. The column titled 'coarse mesh' lists the number of 
triangles in the mesh in Fig. 2.2 that have quality in the range specified in the first 
column. The quality ranges from 2.0 to 5.8. The column titled 'refined mesh' lists 
corresponding values for the mesh determined using a self-similar refinement of the 
background mesh shown in Fig. 2.1a. In this case, the quality ranges from 2.0 to 8.8. 



Range of metric 


coarse mesh 


refined mesh 


2.0-2.4 


1530 


3441 


2.4-2.8 


45 


75 


2.8-3.2 


16 


42 


3.2-3.6 


9 


22 


3.6-4.0 


6 


11 


4.0 - 4.4 


7 


7 


4.4-4.8 


5 


10 


4.8-5.2 


4 


3 


5.2-5.6 


2 


2 


5.6 - 6.0 


2 


3 


6.0 - 6.4 





2 


6.4-6.8 





4 


> 6.8 
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small features or undergoes topological changes, it may no longer belong to V{Th)- 

We illustrate this idea with the example shown in Fig. 2.3. The domain is a 
three-blade propeller that rotates about an axis perpendicular to its plane and passing 
through the center. The background mesh shown in Fig. 2.3a consists of only acute 
angled triangles. It is refined at the center and along the tips of the blades to resolve 
the larger curvatures there. This mesh yields a conforming mesh for every orientation 
of the propeller; a few are shown in Fig. 2.3b. This background mesh is hence a 
universal mesh for each configuration of the propeller. 

An important question in practice is knowing when Q, belongs to V(Th)- Precisely 
characterizing V{Th) is presumably very difficult. Fortunately, it is not essential. 
Rather, the key step in checking if il G T>{Th) is knowing if Th is sufficiently refined in 
the vicinity of d£l. For this, we will require good, computable and local estimates for 
the mesh size in order for the meshing algorithm to succeed. As a step in this direction, 
in [21] we provided upper bounds for the mesh size to guarantee a parameterization 
of dVL over positive edges in Th ■ An analysis of the meshing algorithm is required to 
derive a similar bound for the required mesh size of the background mesh. 

2.4. Details for the implementation. An implementation of the meshing 
algorithm is provided in appendix A. We discuss a few details here, 
(i) Identifying vertices in f2: The first step in Table 2.1 requires identifying 
which vertices of Th lie in f2. This is simplest when Q is represented implicitly, 
as f2 = {x e W 2 : ^(x) < 0}. For then, a vertex v in Th belongs to iff 
^(v) < 0. If such a level set function ip is not known a priori, it can be chosen 
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to be the signed distance function to dfl, (j) :R 2 defined as 



4>(x) = 



— distance(x, 9f2) if x £ Cl, 
distance^, otherwise. 



(2.4) 



Closest point projection: Mapping vertices of positive edges onto 9f2 requires 
computing the closest point projection to <9f2, ir : R 2 — > dfl defined as 



For C 2 -regular domains, ir is related to the signed distance function <fi by 



sufficiently close to <9f2, see [21, Theorem 2.2]. Observe that in the meshing 
algorithm, ir needs to be evaluated only over positive edges and that these edges 
are by definition within a distance h from dil. Hence relation (2.6) can be used 
to compute n if Th is sufficiently refined. This in turn requires computing <f> 
and its derivatives close to the boundary. We refer to appendix A in [23] for a 
discussion on computing (f>, ir, and their derivatives for parametric and implicit 
representations of dfl. 

Relaxing vertices away from dfl: In step 4 in Table 2.1, vertices in f2 that lie 
close to dfl are perturbed away from the boundary. While such perturbations 
can be realized in numerous ways, we have adopted the map in (2.1). Close 
to dtt, V(j)(x) equals the unit outward normal to dfl at n(x). Hence ph(%) 
indeed perturbs vertices away from the boundary. By selecting r to be 0(h) 
in the definition of phi only a small number of vertices near the boundary are 
perturbed. Such a scaling is also essential in the definition of p h because V</> 
may be defined only in a small neighborhood of dSl. In our examples, we pick 
■q ~ 0.3 and r ~ 3ft. Observe that ph docs not move vertices that lie on dCl. 
Hence steps 3 and 4 in Table 2.1 move exclusive sets of vertices. These two steps 
can therefore be performed in cither order. 

Tolerances and Round-off: In identifying which vertices lie in il (step 1 in 
Table 2.1), the effect of tolerances and round-off errors is perhaps unavoidable. 
As a result, a vertex in Q may be (mis) identified as lying in O c and vice versa. 
The effect of tolerances can in fact be understood as introducing small pertur- 
bations in the boundary. Incorrectly identifying vertices in O will change the 
collection of positively cut triangles, positive edges and hence the resulting mesh 
for f2. However, the resulting mesh will be valid provided conditioning angles 
remain acute. In particular, if triangles in the vicinity of <9f2 are acute angled, 
the choice of tolerances and the effect of round-off errors is not critical. The 
resulting mesh may depend on their choice but will be valid nonetheless. 
Background meshes: For a given curved domain fi, assumptions (b)-(d) in 
§2.1 impose restrictions on the background mesh Th- Since the polygon triangu- 
lated by Th is quite arbitrary, assumption (b) is easily satisfied. A simple way to 
satisfy the acute conditioning angle requirement (c) is to ensure that triangles in 
the vicinity of dtt are acute angled. An even simpler way is to use a background 
mesh of all acute angled triangles. For instance, use a mesh of all equilateral 
triangles as done in the example in Fig. 2.1. 

In practice, it is desirable to use an adaptively refined background mesh, de- 
pending on the geometric features of the boundary or on the solution being 



w(x) = arg min distance(x, y). 



(2.5) 



7r(x) = x — (f>(x)V<p(x) 



(2.6) 
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(a) (b) 

Fig. 3.1: Mappings from triangles in a background mesh to curved triangles that 
exactly conform to an immersed domain. As shown in (a), the mapping y>\ takes 
positively cut triangles to curved ones that conform to the boundary exactly. This is 
achieved by ensuring that the restriction of <p\ to a positive edge equals the closest 
point projection 7r. Triangles with two or more vertices in O, i.e., in T^' 1 , are mapped 
affinely. Figure (b) shows part of the curvilinear mesh obtained by using (p 1 ^ in the 
example in Fig. 2.1a. 



approximated. A convenient way of doing so is by triangulating adaptively re- 
fined quadtrees. Bern et al. [4] provide stencils of acute angled triangles to 
tile quadtrees. The interior angles of triangles in these stencils lie between 36° 
and 80°. Therefore, the resulting background meshes automatically satisfy the 
acute conditioning angle requirement. The background mesh shown in Fig. 2.3a 
was constructed in this way. We refer to [23] for more examples of background 
meshes constructed from adaptively refined quadtrees. 

3. Exactly conforming curved elements. Curvilinear discretizations provide 
high-order accurate approximations for curved domains, compared to polygonal ones 
that results from the meshing algorithm discussed above. Curved finite elements 
constructed using such discretizations are indispensable for optimal accuracy with 
high-order interpolations. Curvilinear discretizations broadly fall into two categories. 
In the first kind, curved triangles conform exactly to the domain. In the other, 
curved triangles approximate the domain sufficiently well and are usually defined via 
isoparametric mappings. We consider the former here and the latter in §4. 

Constructing mappings from straight to curved triangles, even with a conforming 
mesh, is a delicate task because there are two conflicting requirements. The resulting 
curved triangle should approximate the domain well. Yet, it should be a sufficiently 
small perturbation of the straight one if interpolation estimates on the latter are 
expected to translate into optimal ones over the curved triangle, see [9]. Below we give 
one such mapping, which generalizes the ones in [13, 28] to the case of nonconforming 
background meshes. 
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3.1. Exactly conforming triangles. Defining a curvilinear mesh that con- 
forms exactly to fi requires only a subtle modification of the meshing algorithm — in- 
stead of mapping vertices of positive edges onto <9f2, we map positive edges themselves 
onto dil. Analogous to the mapping M\ in (2.3) that defined the meshing algorithm, 
we construct a mapping tp 1 ^ triangle- wise over the collection 7^ 0,1 ' 2 . To this end, con- 
sider K G T^ 0,1 ' 2 with vertices {u,v,u>} ordered such that (j){u) > <p{v) > (j){w). For 
K G 7^ 0,1 , set p>\ := M^. Over positively cut triangles K G Tfi , define 

<Pk( x ) '■= /-i 1 , x [XvK (A„ u + (1 - A„) v) + X u X w n(u)] 

[A„7T ((1 - A„) u + X v v) + X v X w tt(v)] 



2(1 -A„) 

+ X w p h (w). (3.1) 

Note that as in (2.3), the dependence on x in (3.1) is implicit in the barycentric 
coordinates X U ,X V and A^,. Unlike in (2.3) however, ip h K in (3.1) is no longer 
afiine over positively cut triangles. Fig. 3.1a depicts the action of <p>\ on triangles in 
7^ 1,2 . Fig. 3.1b shows part of the curvilinear mesh obtained by using the map p>\ 
in the example in Fig. 2.1. Since <p\ differs from only over triangles in 7^ 2 , the 
curvilinear mesh in Fig. 3.1b differs from the mesh in Fig. 2.2b only over the 234 
positively cut triangles. 

Let us examine the definition of p\ for K G in (3.1). By the assumed ordering 
of vertices, the edge uv joining vertices u and v is the positive edge of K. On this 
edge A„, = and X u + X v = 1. So 

ip h K {x Guv) = ^7r((l - X v )u + X v v) + ^n(X u u + (1 - A„) v) = n(x). (3.2) 

Hence ip 1 ^ equals the closest point projection over the positive edge uv. This shows 
that (fix maps the positive edge onto d£l, as depicted in Fig. 3.1a. On the edge ww, 
X v = and A„ + A^, = 1. Then (3.1) reduces to 

ip h K {x G ww) = ^^~Y^ n ( u ) + \ n ( u ) + X w ph(w) = X u tt{u) + X w p h (w), (3.3) 

which is an affine map. Similarly, 

Lp h K (x G WW) = X v ir(v) + X w p h (w). (3.4) 

Eqs.(3.2), (3.3) and (3.4) show that p\ can be interpreted as the interpolation to K, 
of a map that equals 7r on the positive edge and is affine on the remaining two. This 
point of view is also adopted in [13, 17]. For this reason, mappings such as tp\ in 
(3.1) are also commonly termed blending maps and transfinite interpolations. 

Remark: We ought to mention an alternate construction for mapping positively 
cut triangles to curved ones that explicitly uses the meshing algorithm as an interme- 
diate step. In such a construction, the domain O is first meshed using the algorithm 
in §2 and the resulting mesh is then transformed to a curvilinear one that conforms 
to tt. More precisely, with K s := M%(K), the mapping tp^ defined as 

^•~U S °M£ if^GT, 2 (3 ' 5) 
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Fig. 3.2: Defining high-order finite elements by constructing mappings from the ref- 
erence element K to the curved element K c . The figure shows two such mappings. 
The mapping on the left, namely <p\ o A K , first transforms K to the positively cut 
triangle K in the background mesh and then uses <~p\ to map K to the curved triangle 
K c . The second construction on the right, given by p>\ s o Mk ° Ak (also equal to 
^oijf), uses the meshing algorithm as an intermediate step. The reference triangle 
is mapped to K with Ak-, then K mapped to the triangle Kg = Mk(K) which is 
finally transformed to K c with the map p>\ s . 



maps triangles in the collection T^' 1 ' 2 to a curvilinear mesh that conforms exactly 
to 0. Of course, the maps ip>\ and 4>k differ only for positively cut triangles. For 
if £ T ft 2 , observe 

(i) that ipx first maps if to a conforming triangle K$ and then transforms K$ to 
a curved triangle, and 

(ii) that even though (fi h K {K) = ^(K) (as sets in W 2 ), ip h K ^ ip^ in general. The 
two will however be close in a pointwise sense. 

The distinction between the <p^- and V'x for positively cut triangles is illustrated in 
Fig. 3.2. 

The conditions in §2.1 suffice for the mappings ip 1 ^ and to be well defined and 
bijective. The mesh size required near the boundary may be smaller for curvilinear 
discretizations of O using ip 1 ^ or -0^ compared to the mesh size required for the 
meshing algorithm. This is because we require the Jacobian of these maps to be 
positive over the entire element. 

3.2. High-order finite elements. The mappings p h K and from straight 
to curvilinear triangles facilitate a natural construction of curved Lagrange finite 
elements. The idea behind the construction is illustrated in Fig. 3.2. 
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3.2.1. Curved finite elements. Introduce the reference triangle K C R 2 and 
the finite element triplet (K,N k ,V k ). As usual, V k is the set of polynomials over K 
of degree at most k and N k — {N a } a is the set of shape functions that constitute 
a basis for P fc . Associated with N k are the nodes {z a } a € K which are such that 

N a (zb) = S ab . 

Let Ak ■ K — ► K be an affine map from K to K. The curved finite element 
corresponding to K e T° X2 derived from (K,N k ,V k ) is denoted by {K c ,N k ,P k ), 
where K c = tp h K {K) = tp h K {A K {k)) and 

P k = {poA K 1 o( i p h K )- 1 : pe¥ k }. (3.6) 

In particular, shape functions {N a } a over K c are defined by the relation N a o ip 1 ^. o 
Ak = N a . Nodes {z^} a in the curved element are located at z c a = ip K (AK(z a ))- 
Hereafter, the set {i a } will be chosen so that finite element functions over ft are in 
C°. 

The choice of ip K over tp K to define the curved element (K c ,N k ,P k ) above was 
arbitrary; replacing each instance of (p K by tp K yields a curved element as well. In 
fact, it may be more convenient to incorporate tp K into existing finite element codes 
based on conforming meshes. In the following example as well as the ones in §3.3, we 
have adopted the curved elements based on the map ip K . 

3.2.2. Optimal convergence: numerical example. We demonstrate opti- 
mal convergence using the curved finite elements described above with a numerical 
example. Although the given construction for curved elements is a standard one, 
the example helps show that the mapping ip K in (3.1) satisfies the conditions in [9] 
for optimal interpolation estimates over the curved element. We consider the model 
problem 

Au = in Q = {r = ^/x 2 + y 2 < 1}, (3.7a) 
u = e y sin x on dtt. (3.7b) 

The solution to (3.7) is the smooth function u(x,y) = e v smx. The weak form of 
(3.7) is to find u e = {v e H 1 ^) : v\ gn = e^sinrr} such that 

/ Vu-Vvdn = VveH^fl). (3.8) 
Jn 

To compute finite element approximations u h of u, il is immersed in background 
meshes of equilateral triangles. The coarsest background mesh has mesh size ft-o — 
0.27. Fig. 3.3 shows the convergence of the solution computed with standard La- 
grange elements (over K), as the background mesh is refined (/i-refinement). Dirichlet 
boundary conditions were imposed by interpolating the prescribed function in (3.7b) 
at the nodes of curved elements lying on the boundary. We used sufficiently accurate 
quadrature rules to evaluate the stiffness matrix, see §4.1. The convergence rate in the 
L 2 (f2)-norm is optimal for linear, quadratic, cubic and quartic elements (k = 1,2,3 
and 4 respectively). 

Examining ||u — Uh\\L 2 (n) m Fig- 3.3 also reveals that for a given background 
mesh, the error decreases with the element order k. This demonstrates that the 
curved elements are well suited for p-refinemcnt — progressively accurate solutions can 
be computed by just increasing the element order while using the same background 
mesh. Moreover, the data in Fig. 3.3 shows that the error is 0(h k+1 ), which is 
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Fig. 3.3: Optimal conver- 
gence of the finite element 
solution Uh, computed using 
exactly conforming curved 
elements, to the exact one 
u of problem (3.7). The 
plot shows convergence in 
the L 2 (f2)-norm as the back- 
ground mesh is refined. The 
rate of convergence is opti- 
mal for linear, quadratic, cu- 
bic and quartic elements. 



optimal in the element order k for each given (sufficiently small) mesh size h of the 
background mesh. Such optimal convergence rates would also be obtained with the 
isoparametric curved elements described subsequently in §4. 

3.3. Applications: Evolving fluid domains. Next, we present two applica- 
tions using the curved elements described above. In both examples, a fixed background 
mesh serves as the universal mesh for an evolving fluid domain. 

3.3.1. Flow with a rotating component. We consider the example mentioned 
in §1, of a propeller mixing fluid in a closed container. The problem setup is the same 
one illustrated in Fig. 2.3, although the container B is larger. The propeller P is 
assumed to be rigid and impermeable. It rotates with constant angular velocity lo 
about an axis passing through its center and perpendicular to its plane. The fluid in 
the container is incompressible and has viscosity fi. Its kinematics is governed by the 
familiar equations for Stokes flow, 

^div(Vu t ) = Vpt, (3.9a) 
div (u t ) = 0, (3.9b) 

relating the flow velocity u t and pressure p t at time t. No-slip boundary conditions 
along the walls of the container and the boundary of the propeller imply 

JO on as, 
u t = < (3.10) 
[ruje g ondP t , 

where P t is the configuration of the propeller at time t and {e r ,eg} is a canonical 
polar basis for a polar coordinate system with origin at the center of P. Since (3.9) 
and (3.10) determine the pressure only up to a constant, pt is assigned to be zero at 
one point in the flow domain, i.e., we set p t = at x n e B \ P t . 

A triangulation of B, similar to the one shown in Fig. 2.3a, serves as the universal 
mesh for the fluid domain B \ P t . The mesh is refined further near the tips of the 
blades to resolve features of the flow there. We adopt (curved) Taylor-Hood elements 
for the finite element solution of (3.9), i.e., the element (K C ,N 2 ,P 2 ) for the velocity 
u t and (K c , N , P 1 ) for the pressure p t - See [11] for a discussion on the Taylor- 
Hood clement and for the weak form of this problem, which we have omitted here. 
Dirichlet boundary conditions are imposed by interpolating (3.10) at the nodes lying 
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(b) 




(c) (d) 

Fig. 3.4: Simulating Stokes flow driven by a rotating propeller in a closed container. 
At each instant, the same background mesh is used to determine a curvilinear mesh 
that conforms exactly to the fluid domain. Figures (a) and (b) show the elements 
close to the propeller in the resulting mesh at two distinct times. Contours of the 
horizontal component of the flow velocity at these times are shown in figures (c) and 
(d). To highlight the flow pattern, the contours are shown with zebra shading in (d). 



on the boundary. Figs. 3.4a and 3.4b show the curvilinear mesh conforming to the 
fluid domain at three different time instants. Since the mesh is quite refined, only 
the elements near the propeller are shown. Corresponding to these orientations of 
the propeller, Figs. 3.4c and 3.4d show contours of the horizontal component of the 
velocity computed with fi = 0.01 and lj = 2. 

3.3.2. Flow interaction with a rigid disc. In the second example, we consider 
the interaction between a fluid and a rigid solid. The problem is to determine the 
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Fig. 3.5: Figure (a) shows the initial setup for the problem of a rigid disc D interact- 
ing with an incompressible fluid. The disc is attached to the origin by a linear spring. 
Inflow and outflow boundary conditions for the flow through the channel B are indi- 
cated. A simple, unstructured mesh of B serves as the universal mesh for the fluid for 
the entire duration of the simulation. Figure (b) shows the trajectory computed for 
the disc as it moves from its initial position to an (approximate) equilibrium position. 



trajectory of a rigid disc D of radius R and mass m immersed in an incompressible, 
viscous fluid flowing through a square shaped channel B of side L. The disc is attached 
to the origin O located at the mid-point of the left end of the channel by a linear spring 
with spring constant k and equilibrium length £ . The problem setup is shown in Fig. 
3.5a. 

We assume that the kinematics of the fluid is governed by the equations for Stokes 
flow given in (3.9), retaining the notation introduced there. In a Cartesian coordinate 
system centered at O, inflow and outflow boundary conditions are prescribed at the 
two ends of the flow channel as 

Ut = (L-2|y|)e x if* = 0,1,. (3.11) 

Denote the position of the center of the disc at time t by c(t) and the disc centered 
at c(t) by D t . No-slip boundary conditions along the horizontal walls of the channel 
and along the boundary of the disc imply 

fo if|y|=L/2, 
I c(i) on oD t . 

Force balance for the disc is given by 

mc = fc(l — ^r]c+ f (Tf-n t ds, (3.13) 

V \ C \J JdD t 

where n t is the unit outward normal to dD t and the stress 07 in the fluid is computed 



as 

07 = -p t I + A* (Vu ( + Vu 



17 



Fig. 3.6: Contours of 
the horizontal com- 
ponent of the flow 
velocity at a non- 
equilibrium position 
of the disc in the 
problem described in 
Fig. 3.5a. 

Balance equations (3.9), (3.13), boundary conditions (3.11), (3.12), initial con- 
ditions c(0) = Co, 6(0) = 0, and pt(xo € B \ D t ) = together constitute a coupled 
system of equations for the unknowns (u t ,p t ) and c(t). We use (curved) Taylor-Hood 
elements for the flow solution as before and adopt a staggered time integration scheme. 
At each time instant t n , given the center of the disc c„ and its velocity c„, we define 
curvilinear elements for the flow variables over B \ D t . We then compute the flow 
solution (ujj,p^) at this time. The net force on the disc is evaluated using (3.13). 
Using central differences, we update the position and velocity of the disc to the next 
instant t n+ \ and repeat the process. 

The background mesh shown in Fig. 3.5a serves as a universal mesh for the flow 
domain B\D t . Fig. 3.5b shows the trajectory determined for the disc and its final 
configuration computed with parameters L — 1 for the container, /i = 0.01 for the 
fluid, R = L/10,to = 1 for the disc, k = l,£o = L/A for the spring, c = (0,L/4) for 
the initial position of the disc and a time step At = 0.05. At an intermediate position 
of the disc, contours of the horizontal component of the flow velocity along with a few 
stream lines are shown in Fig. 3.6. The trajectory of the disc plotted in Fig. 3.5b 
shows that the disc eventually settles to an equilibrium position balancing the forces 
exerted by the fluid and the spring. 

We conclude this section mentioning that there are various numerical methods in 
the literature for such problems over changing flow domains, see for instance [3, 10, 
24, 27]. 

4. Isoparametric mappings. Isoparametric mappings provide a convenient 
way of approximating curved domains with a desired accuracy. A systematic defi- 
nition of these polynomial mappings results naturally from interpolating Lp\ o Ak or 
*Pk ° A K at selected points. Since tp^ o A K and tp 1 ^ o A K are affine for K e T°'\ 
these maps differ from their interpolants only for positively cut triangles (7^ 2 ). With 
an isoparametric map, positive edges are mapped to curved ones that interpolate the 
boundary at a few points. Fig. 4.1 depicts this for the case of a quadratic element. 

Following the notation in §3.2, introduce the interpolation operator II fe : / g 
[C°(K)] 2 — > J2 a f(za)N a £ P k . The isoparametric map over K corresponding to 
triangle K <G T^' 1 ' 2 is defined as 

I h K :=ti k ( V h K oA K )^Y.^ h K^) N - C 4 - 1 ) 

a 
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(a) 



Z3 
(b) 



z 3 
(c) 



Fig. 4.1: Isoparametric mappings for positively cut triangles in the background mesh. 
The reference quadratic clement is shown in (a) . The isoparametric maps constructed 
by interpolating the mappings ip\ o Ak and tp K o Ak are shown in (a) and (b). The 
curved boundary of the element interpolates the boundary at the points where the 
nodes z\, z% and £4 are mapped. Notice that the mid-side node £4 is generally mapped 
differently by ip K and tpK- 



Using xp K instead of ip\ yields a different isoparametric map 

J h K := Tl k (iP h K oA K )=Y J V h Ks {M K {z a ))N a . (4.2) 

a 

As illustrated in Fig. 4.1 for the case of a quadratic element, I K 7^ J K in general. In 
the figure for instance, I K (zi) = n i( z i + z 2)/2) while J K {z4) = it((it(zi) + n(z2))/2). 
Nonetheless, the two maps will be close for small values of h. Curved finite elements 
using these isoparametric maps are constructed just as in §3.2 by replacing maps tp K 
or i\) K by their respective interpolants. 

Compared to the exactly conforming curved elements, isoparametric elements 
require fewer evaluations of ir in general. Notice from (3.1) that once w is computed 
at the vertices of the positive edge, defining I K (or J K ) requires evaluating 7r at most 
twice per node in the clement. In contrast, computing <p\ (or ip^) at each quadrature 
point requires two new evaluations in the conforming curved element. Furthermore, 
computing derivatives of shape functions in the isoparametric element does not require 
computing derivatives of n. However in the conforming curved element, shape function 
derivatives depend on Vf 1 ^ (or V?/^) which in turn depend on Vn. 

4.1. Quadrature for curved elements. For optimal convergence and accu- 
racy of numerical solutions computed using curved elements, we naturally require 
sufficiently accurate quadrature rules for integration over curvilinear domains. Fol- 
lowing standard practice, these quadrature rules need only be defined over the ref- 
erence triangle, since integrals over a curved element K c can be performed over K 
using the correspondence provided by the mappings l\ (or J^-). We adopted for 
curved elements the same quadrature rules needed for straight elements, as explained 
in [8]. For example, a quadrature rule that exactly integrates quadratic polynomi- 
als over the reference element suffices for isoparametric quadratic elements, so three 
quadrature points were adopted for those. We have used such integration rules in all 
of our examples, including the ones with exactly conforming curved elements. These 
examples suggest that the quadrature rules for straight elements are also enough to 
obtain optimal convergence rates with exactly curved elements. We also note that 
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different rules can be used for integrating each term in the weak form of a problem, 
for instance the mass and stiffness matrices, and the force vector. 

4.2. Circular plate in bending. We consider the problem of a thick, circular, 
clastic plate bending under the action of a uniform external load. As pointed out in [2] , 
this problem highlights the importance of representing curved boundaries accurately. 
We present this example to emphasize the distinction between representing a curved 
boundary exactly using the elements in §3 and approximating it with isoparametric 
elements as defined above. Consider a Cartesian coordinate system (x, y, z) with basis 
{e^, e y , e z }. The domain of the problem is the set £1 x (—t/2, t/2), where fHs a circle 
of radius R = 3.142 centered at the origin and contained in the plane of e x , e y , and t 
is the thickness along e z . We assume that the displacement field u of the plate is of 
the form 

u(x, y, z) = -z (i9 x (x, y)e x + tf y (x, y)e y ) + w(x, y)e z , (4.3) 

which corresponds to the Reissner-Mindlin model for a thick plate in bending, see 
[2, 7]. In (4.3), w is the transverse displacement of points in the mid-plane il while fix 
and fly represent the infinitesimal rotations of fibers normal to the mid-plane about 
the axes e y and e x , respectively. We consider a "soft-simple support" for the plate, 
which implies the boundary conditions 

w = Ol 
0-t = 0] 

where t is the unit tangent to <9f2 and i? = ("& x , "& y ). The plate is loaded by a constant 
force 2p normal to its top face f2 x {t/2}. The elasticity problem is then to find 

u e + we z : t? e Hj.toe flo( fi )} - 

where Hj(fi) := {■& G [H 1 ^)} 2 : ■& ■ t = on dtt}, 

which minimizes the strain energy functional 

J[u] = \ f {A[tr(e(u))] 2 + 2 M£ (u) : e(u)} - / pw, (4.5) 

where A, \x are material parameters called Lame constants and e(u) = (Vu + Vu T )/2 
is the usual infinitesimal strain tensor. Introducing assumption (4.3) in (4.5) and 
integrating along the thickness reduces (4.5) to a problem over f2: find ( l d,w) € 
Hj(f2) x -ffo(Sl) that minimizes the functional 

F[(tf, w)} = l -J^ {A [tr (e(m 2 + ■ 

+ ^/ n ll^-V^|| 2 -f J a PW. (4.6) 

We compare the transverse displacements at the center of the plate as the back- 
ground mesh for is refined, while using curved quadratic elements. We pick A = 
H = l,t = R/A for the plate and p = lxlCP 3 for the loading. We take as a refer- 
ence value wo — 5.3407075 xl0~ 2 , computed with exactly conforming curved quartic 
elements (see §3.2) and a refined background mesh of equilateral triangles. Table 
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(4.4) 



Tabic 4.1: Transverse displacements at the center of the circular plate 0, computed 
with curved quadratic elements. The reference value is wq = 5.3407075x 10~ 2 . The 
columns titled 'conforming' and 'isoparametric' list the values computed with exactly 
conforming elements and isoparametric elements, respectively. The column 'modified 
isoparametric' contains the values computed with isoparametric elements but while 
imposing the constraint "d ■ t = using one of the two possible tangents t at vertices 
on the curved edges of the element. The coarsest mesh size of the background mesh 
is h ~ 0.28i?. 



mesh size 


conforming (xlO 2 ) 


isoparametric (xlO a ) 


modified isoparametric (xlO a ) 


h 


5.3370662 


2.2343001 


5.3383666 


ho/2 


5.3401468 


1.9557577 


5.3400195 


hoi '4 


5.3406606 


1.8019916 


5.3406438 


h /8 


5.3407108 







4.1 lists the displacements computed with exactly conforming quadratic elements and 
with isoparametric quadratic elements. 

From the table, we see that the displacements computed with the conforming 
elements converge to wq. But somewhat surprisingly, those computed with quadratic 
isoparametric elements fail to even come close to w - This is a consequence of enforcing 
the constraint on rotations in (4.4) on the approximate curved boundary realized with 
isoparametric elements. The unit tangent to this boundary fails to be continuous at 
vertices that lie on it. Consequently, rotations equal zero at each vertex on the 
boundary. With exactly conforming curved elements, 9f2 is represented exactly and 
this issue is avoided. 

The above discussion shows that the constraint in (4.4) needs to be enforced 
differently. For instance, we could select one of the tangents at each vertex on the 
approximate boundary to enforce the constraint i? • t = 0. The displacements com- 
puted with quadratic isoparametric elements by enforcing the constraint in this way 
are listed under the column title 'modified isoparametric' in Table 4.1. These values 
are clearly more accurate and converge to wo- 

Remark: We have used the same finite element spaces for both transverse displace- 
ments (w) and rotations (d x ,-d y ) in the above calculations. It is well known that for 
thin plates (t <C R), such a choice of spaces in the Reissner-Mindlin model results in 
locking, cf. [7] . To avoid adopting very small mesh sizes for accuracy, we deliberately 
chose a large thickness t = R/4 in the example. 

5. Rationale behind algorithm. The meshing algorithm determines a con- 
forming mesh for by perturbing vertices in the triangulation 7^ ' . Here we briefly 
discuss the rationale behind the vertex adjustments we perform — why we perturb 
vertices in a particular way and when such perturbations are possible. 
Mapping vertices to their closest point on 9fi: Each vertex of a positive edge 
is mapped to its closest point on the boundary. This step transforms T^' 1 ' 2 into a 
conforming mesh for fL The reason is an intuitive one. As we discuss below, the only 
vertices in T^' 1 ' 2 that lie outside f2 are the vertices of positive edges. Moreover, the 
collection of positive edges are the boundary edges of T^' 1,2 because each of these 
edges belongs to just one positively cut triangle. By snapping vertices of positive 
edges onto dil, every vertex in the resulting mesh belongs to Q and the boundary 
edges in the final mesh interpolate dfl. 
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Fig. 5.1: Vertices in near its boundary need to be relaxed away from the boundary in 
the meshing algorithm to ensure good quality triangles. The figure shows a positively 
cut triangle with vertices {u, v,w} and positive edge uv. Step 3 in Table 2.1 moves 
vertices u and v to n(u) and tt(v) respectively. The quality of the resulting triangle 
with vertices {tt(u), tt(v), w}(shaded in gray) depends on the distance of w from <9f2. 
The three cases in the figure show that the closer w is to <9f2, the poorer the quality 
of the shaded triangle. Perturbing w away from dfl alleviates this. 



A small mesh size near dfl is essential for this perturbation step. For otherwise, 
7T may be multi-valued at vertices of positive edges. A small mesh size is also required 
to know that positive edges are boundary edges of 7^°' 1,2 . If the mesh size is too 
large, it is possible for two positively cut triangles to share a common positive edge. 
The acute conditioning angle assumption in §2.1 is critical as well — it ensures that 
both vertices of a positive edge never map to the same point in <9f2 hence preventing 
positively cut triangles from being mapped to degenerate ones. In fact, using the 
assumptions in §2.1, we proved in [21] that the restriction of tt to the collection of 
positive edges is a homeomorphism onto dSl with Jacobian bounded away from zero. 
Hence we know that moving vertices of positive edges onto the boundary does not 
yield a tangled mesh. The bound for the Jacobian also implies that positive edges are 
mapped to interpolating edges whose lengths are neither too small nor too large. We 
refer to [21, 23] for a detailed discussion on the sufficiency of the assumptions in §2.1 
to show that <9f2 is parameterized over the collection of positive edges. There we also 
mention a way of relaxing the acute conditioning angle assumption. 

Relaxing vertices away from dtt: Next we explain using an example why it is 
necessary to relax vertices away from the boundary (step 4 in Table 2.1). Consider a 
positively cut triangle K with vertices {u, v, w} and positive edge uv as shown in Fig. 
5.1. By snapping vertices u and v onto d£l in step 3 in the algorithm, K is mapped 
to a triangle K with vertices {tt(u),tt(v),w} (shaded in gray in the figure). We know 
from the above discussion that the length of the edge tt(u)it(v) cannot be too small. 
However, lengths edges ir(u)w and ir(y)w can be arbitrarily small. As depicted in the 
figure, the closer w is to dCl, the poorer the quality of triangle K. To alleviate this, 
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we move vertex w away from dfl by a small distance. In turn, to accommodate such 
perturbations, vertices in a small neighborhood of dfl are relaxed away from it. 

In the map p h used to relax vertices in the algorithm, the parameter r determines 
the neighborhood of <9f2 in which vertices are perturbed. Since we pick r to be a few 
multiples of the mesh size near dQ, ph is only a local perturbation near dfl. To prevent 
this step from inducing overlapping or poorly shaped triangles, a small-enough mesh 
size is also needed. For example, it needs to be well defined in the r-neighborhood 
of dQ, so that ph is well defined as well over the vertices that are relaxed. Although 
a detailed analysis is still needed, this step does not appear to pose more stringent 
requirements on the mesh size that those posed by y>\ or 

6. Concluding remarks. The meshing algorithm, the mappings to curvilinear 
triangles used in constructing high-order curved finite elements, and the idea of uni- 
versal meshes are useful tools for an important class of computationally challenging 
problems. By employing them in problems with evolving fluid domains while using a 
single background mesh, we demonstrated the algorithmic advantages they offer. We 
envision their application to more demanding problems ranging from dynamic crack 
propagation to phase transformations. 

For these tools to also be useful in realistic engineering applications, important 
questions remain. A significant one is knowing what is a sufficiently small mesh size 
for the background mesh. A computable estimate is valuable because it can help 
determine if and when a background(univcrsal) mesh needs to be changed during 
the course of simulating an evolving domain. The mesh size estimates in [21] for 
parameterizing the immersed boundary will be useful in determining such bounds. 

A second challenge is summarized by the fact that while we can mesh extremely 
complex smooth domains using a simple background mesh, we have not specified how 
to mesh a square. This is essentially a consequence of choosing the closest point pro- 
jection to parameterize the boundary. Whether domains with corners, cracks, and 
interfaces can be handled without introducing additional restrictions on the back- 
ground mesh remains to be seen. An important step towards meshing such domains 
is parameterizing immersed curves with end points and corners. We have shown how 
to do this in [23]. 

We think the ideas introduced here can be extended to meshing three dimensional 
domains immersed in background meshes of tetrahedral meshes. An analysis will 
reveal the necessary requirements on the background mesh. 

Finally, we mention that maintaining the regularity of evolving boundaries has 
been a recurring challenge in numerical methods for moving boundary problems. It 
requires careful choices for representing the domain and for schemes to advance the 
boundary. The literature on these topics is growing and will continue to benefit from 
each new contribution. 
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Appendix A. Implementation for the meshing algorithm. We provide 
a simple implementation of the meshing algorithm in §2 to determine a conforming 
mesh for Q,. We assume that the background mesh Th is specified by (i) a list of 
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coordinates V of its vertices, (ii) a numbering 7 for the vertices, and (iii) a list of 
triangle connectivities C which are 3-tuples of vertex numbers. Vertex with number 
i g 7 is denoted Vi G V. A triangle with vertices {vi, Vj,Vk} has connectivity (i, j, k) g 
C. The algorithm returns a conforming triangulation for Q specified by a set of vertices 
V n , a numbering I n C 7 for these vertices and a connectivity list C n C C for triangles 
in the mesh. 

Require: vertex coordinates V, vertex numbering 7, triangle connectivities C. 
Require: Choose r\ g (0, 1) and r equal to a few multiples of h 

l: Initialize 7° <- 0, <- 0, C n <- 

2: for alii g I do 

3: if Vi g n then 

4: Si <— true 

5: else 

6: Si 4— false 

7: end if 

8: if Si then 

9: Append i to I Q 

10: if — r < <j){vi) then 

11: Append p/ l (w J ) to V n 

12: else 

13: Append vi to F° 

14: end if 

15: end if 
16: end for 

17: 7+^0 

18: for all k) g C do 

19: i_^0,i+^0 

20: for £ g {i, j, fc} do 

21: if se then 

22: Append £ to i_ 

23: else 

24: Append I to i + 

25: end if 

26: end for 

27: if #i_ > 1 then 

28: Append (i, j, k) to C° 

29: if #i + = 2 then 

30: Ensure conditioning angle of triangle (i, j, k) is acute. 

31: for £ g i + do 

32: Append ^ to 7+ 

33: end for 

34: end if 

35: end if 

36: end for 

37: for all i g 7+ do 

38: Append i to 7 n 

39: Append n(vi) to 

40: end for 

41: return triangulation (F n ,7 n ,C a ) for 
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